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Abstract 

Accurate models of carrier transport are essential for describing the electronic properties of 
semiconductor materials. To the best of our knowledge, the current models following the frame¬ 
work of the Boltzmann transport equation (BTE) either rely heavily on experimental data (i.e., 
semi-empirical), or utilize simplifying assumptions, such as the constant relaxation time approxi¬ 
mation (BTE-cRTA). While these models offer valuable physical insights and accurate calculations 
of transport properties in some cases, they often lack sufficient accuracy - particularly in capturing 
the correct trends with temperature and carrier concentration. We present here a transport model 
for calculating low-field electrical drift mobility and Seebeck coefficient of n-type semiconductors, 
by explicitly considering relevant physical phenomena (i.e. elastic and inelastic scattering mecha¬ 
nisms). We first rewrite expressions for the rates of elastic scattering mechanisms, in terms of ab 
initio properties, such as the band structure, density of states, and polar optical phonon frequency. 
We then solve the linear BTE to obtain the perturbation to the electron distribution resulting 
from the dominant scattering mechanisms - and use this to calculate the overall mobility and 
Seebeck coefficient. Using our model, we accurately calculate electrical transport properties of the 
compound n-type semiconductors, GaAs and InN, over various ranges of temperature and carrier 
concentration. Our fully predictive model provides high accuracy when compared to experimental 
measurements on both GaAs and InN, and vastly outperforms both semi-empirical models and the 
BTE-cRTA. Therefore, we assert that this approach represents a first step towards a fully ab initio 
carrier transport model that is valid in all compound semiconductors. 

PACS numbers: 72.20.-i, 73.61.Ey, 31.15.A-, 71.20.Nr 

Keywords: band transport model, mobility, Seebeck coefficient, electron scattering, ionized impurity, acous¬ 
tic phonon, polar optical phonon, band structure, density of states, III-V semiconductors 
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I. INTRODUCTION 


Accurate models of carrier transport are essential for describing the electronic properties 
of semiconductor materials, which are particularly important for clean energy applications 
ranging from photovoltaics to thermoelectrics to photoelectrocatalysts. There has been an 
increased focus on using compound semiconductors, including those that are degenerately 
and heavily doped, for these applications. To better understand existing materials and 
discover new ones, a fully predictive model that correlates electronic structure to properties 
is essential. Unfortunately, to the best of our knowledge, no model, based on ab initio 
calculations, currently exists to fully capture the elastic and inelastic scattering effects of 
charge carriers; as a result, errors arise when utilizing the current models. While an ab 
initio model will certainly improve our understanding of the carrier transport mechanism(s) 
in existing semiconductors, it can also aid in the search for high-performing materials by 
improving the accuracy of high-throughput computation^ 2 ! 

There currently exist two main categories of models, based on the Boltzmann transport 
equation (BTE), for calculating the conductivity and Seebeck coefficient of semiconductors 
that are governed by band conduction. The first category of BTE-based models are com¬ 
monly known as single parabolic band models, even though the treatment of the conduction 
band may not be explicitly parabolic. These models can be described as ” semi-empirical”, 
since experimentally measured parameters, such as the electron or hole effective mass, band 
gap, dielectric constant and polar optical (PO) phonon frequency, are used in closed-form 
expressions for the various scattering rates. Note that the overall mobility due to elastic 
scattering is calculated by averaging, according to Matthiessen’s rule, the mobilities due to 
each scattering contribution. The main adjustable parameter in these models is the effective 
mass, which can be varied to fit the calculated transport properties to the experimental 
measurements. While such models often impressively capture the changes in properties over 
various ranges of temperature and carrier concentration, they are restricted to the materials 
for which experimental data are available; therefore, the predictability of such models are 
very limited. 

There are numerous examples of models in this categorjPP, such as that by EhrenreiclP, 
who modeled the GaAs band structure and PO-phonon scattering by reviewing the exper¬ 
imental dataP, and that by Sankey et alP, who considered the effects of resonance, ionized 
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impurity, and polar optical phonon scattering in GaAs. In these models, all of the scat¬ 
tering mechanisms are commonly treated using the relaxation time approximation (RTA); 
here, the relaxation time is written as a power law function of energy - thus, the details of 
clastic and inelastic scattering (e.g., PO phonon) captured by the ab initio band structure 
are disregarded. Scattering rates, particularly inelastic ones, have already been shown to not 
follow such power law distributions® so the basic assumptions fail. Even in cases where the 
BTE is explicitly solved for PO phonon and the perturbation to the electronic distribution 
is obtained without the RTA assumption, the results are still heavily dependent on available 
experimental data. As an example, Miller et ah 8 used the latter approach to calculate the 
mobility and Seebeck coefficient of InN samples, which had been grown by molecular beam 
epitaxy (MBE) and plasma assisted MBE so that all exhibited heteroepitaxial growth with 
linear charged dislocations; thus, these dislocations were found to be the limiting scattering 
mechanism. 

The second category of BTE-based models relies on the ab initio band structure of the 
material, rather than specific experimentally measured parameters, but generally utilizes 
the relaxation time approximation to the BTE (BTE-RTA) as a simplification. Restrepo et 
alP calculated the mobility of n-doped silicon at different electron concentrations in BTE- 
RTA and ab initio framework where electron-phonon interactions are treated as elastic with 
the electron distribution unchanged from the equilibrium Fermi-Dirac. On the other hand, 
the constant relaxation time approximation (BTE-cRTA) simplifies the equation even more, 
which enables closed form expressions for both conductivity divided by relaxation time and 
Seebeck coefficient. The advantage of these models is the ability to calculate properties 
of new materials, for which experimental data is unavailable. This type of model works 
well for some materials for which the relaxation time is fairly constant, as evidenced by the 
work of Madsen and Singh 10 . However, inelastic scattering mechanisms change the electron 
energy and directly affect the distribution. Lumping all the elastic and inelastic scattering 
mechanisms into a single constant and assuming an equilibrium Fermi-Dirac distribution 
in BTE-cRTA framework greatly damages the predictive ability of such models; as an ex¬ 
ample, transport properties in some cases are very far from experimental measurements. 
Furthermore, the relaxation time constant is usually determined by fitting the calculated 
conductivity to experimental data. It should be noted that the calculation of this constant 
is not necessary when calculating the Seebeck coefficient. This is due to the simplifying 
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assumptions that the relaxation time is both a constant and direction independent 10 which 
does not always hold. Therefore, BTE-cRTA suffers not only from inaccuracy in predict¬ 
ing the changes of properties with temperature or carrier concentration in many materials, 
but also from lack of pure predictability since it still relies on experimental data for the 
computation of the relaxation time. 

Instead, we propose that accurate calculations of electronic transport properties, within 
the Boltzmann transport framework, are possible by combining relevant treatment of the 
elastic and inelastic scattering mechanisms with ab initio calculations of the electronic and 
phonon band structures. Ultimately, an ab initio theory for carrier transport will need to 
qualitatively and quantitatively predict trends in material properties, such as conductivity 
and Seebeck coefficient, as a function of temperature or carrier concentration. Validation 
of the theory against experimentally measured properties will thus give insight into which 
scattering effects are dominant. 

In this paper, we present a band transport model for calculating low-held electrical drift 
mobility and Seebeck coefficient of n-type semiconductors. We then validate our model by 
calculating the properties of two III-V semiconductors, GaAs and InN, with different carrier 
concentrations over various temperatures, and comparing them to experimental values as 
well as those calculated using the other transport models described above. We choose these 
materials because the ab initio band structure of GaAs is similar to those used in the earlier 
semi-empirical models at it can be reasonably well described with a single band model, 
whereas the ab initio band structure of InN and the limiting scattering mechanisms are 
quite different; thus, these two materials allow us to bracket the range of expected behavior 
of our proposed model. 


II. CARRIER TRANSPORT MODEL 


A. Solution to the Boltzmann Transport Equation 


In order to calculate the mobility and Seebeck coefficient, we solve the Boltzmann trans¬ 
port equation (BTE) using Rode’s iterative method 3 8 11 17 (Appendix A 2) to obtain the 


electron distribution in response to a small driving force (e.g. a small electric held or a 
small temperature gradient). It is important to note that we do not use the relaxation time 
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approximation (RTA) in solving the BTE, so neither a variable nor a constant relaxation 
time appears in this expression. Due to the assumption of a small driving force, we aim to 
calculate only the linear response to the perturbation; thus, the general form of the electron 
distribution remains the at equilibrium Fermi-Dirac distribution. We can then write: 

/ (k) = f 0 [e (k)] + xg (. k ) (1) 


where / is the actual distribution of the electrons, including both elastic and inelastic scat¬ 
tering mechanisms, / 0 is the equilibrium Fermi-Dirac distribution, x is the cosine of the 
angle between the small driving force and k, g (k) is the perturbation to the distribution 
caused by the small driving force and finally k = |k|. For the sake of simplicity, we express 
the conduction band as the average energy of the electrons as a function of distance, k, from 
the conduction band minimum (CBM) which is often at the center of the Brillouin Zone (i.e. 
T point); furthermore, we assume that the small driving force is aligned with k (i.e., x—1). 
Although this is similar in spirit to the isotropic band assumption, we take the anisotropy 
into account by averaging the energy values of the ab initio calculated band structure, e (k), 
as a function of k rather than explicitly including k in every direction. Alternatively, if we 
wish to consider the directional transport properties, we can include the calculated band 
structure only in that specific direction. Here, we will focus on calculating and reporting 
the overall average mobility and Seebeck coefficient. 

Our goal is to calculate the perturbation to the distribution^, g(k). In the reformu¬ 
lated Boltzmann transport equation shown in Equation |2j there are scattering-in, S{ (g), 
and scattering-out, S 0 , terms for inelastic scattering mechanisms. However, these terms 
also depend, in turn, on the electronic distribution as well as elastic scattering rates, u e \. 
Therefore, the BTE must be solved self-consistently to obtain g (k): 


9(k) 


Si bWj - v W (If) - f (f) 

S 0 (k ) + u ei (k ) 


( 2 ) 


where E is the low electric field and v (. k ) is the electron group velocity. The derivation of 
the BTE in the form shown in Equation [ 2 ] can be found in the literature 3 . The inelastic 
scattering mechanism that tends to dominate at room temperature is polar optical (PO) 
phonon scattering, for which we have provided the description of the Si (g) and S Q terms in 
Equations |A9| and |A10| The influence of inelastic scattering mechanisms on g, and therefore 
the overall mobility, are captured through the terms S t (g) and S Q in Equation [2j while elastic 
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scattering mechanisms affect the overall mobility by the term v e i. This term is the sum of all 
elastic scattering rates inside the material, it can be evaluated according to Matthiessen’s 
rule: 

Vei (k) = Vu (. k ) + u pe ( k ) + L> de (k) + Vdis (k) (3) 


where the subscripts el, ii, pe, de, and dis stand for elastic, ionized impurity, piezoelectric, 
deformation potential and dislocation scattering rates, respectively. Therefore, the effect 
of relevant elastic and inelastic scattering mechanisms are taken into account by explicitly 
solving the BTE (Equation [2j to obtain g (k). 

When calculating various properties, several terms in Equation [2] will be set to zero. For 
a Seebeck coefficient, S, calculation, the applied electric driving force, — (Ik)’ * s se ^ 

zero. Only the thermal driving force, v (§f), in Equation^ is taken into consideration when 
calculating the perturbation to the electron distribution 3 . Assuming a uniform electron 
concentration over the space at which a small temperature difference exists, the Seebeck 
coefficient isR 

£F I k 2 f (1 - /) (5J5,) dk 


S =^ 

e 


dT 

dz 


( 4 ) 


k B T Sk*f(l-f)dk 
For a mobility calculation, the applied thermal driving force in Equation [2] is set to zero, 
so that only the contribution of the electric driving force is included. The mobility is: 


l M «0 (lY(i)dk 

* 3 /(DT* 

Note that in Equationpl the free electron density of states, , has been used, which would 
limit its applicability in compound semiconductors. Thus, the replacement of this term by 
its ab initio -calculated counterpart would greatly improve the accuracy of the resulting 
mobility. Furthermore, the scalar group velocity, v (k), is used since the energy is averaged 
as a function of distance from the T point. In general, we use the band structure, density 
of state, electron group velocity, conduction band wavefunction admixture and PO phonon 
frequency in calculating the mobility and Seebeck coefficient. Therefore, all of the required 
inputs to Equation [5] are calculated ab initio , which greatly enhances the predictability of the 
model. In other words, the main difference between our proposed carrier transport model 
and previous semi-empirical models 3 6 8 11 14 18 20 is the use of ab initio parameters instead of 
experimentally measured electron effective mass, band gap, etc., which eliminates the need 
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for theories such as k-p to describe the nonparabolicity or anisotropy of the conduction band. 
Instead, for calculating the overlap integral, we express the conduction band wavefunction 
as a linear combination of s-type and p-type basis functions, with coefficients of a and c, 
respectively 3 . These coefficients can be directly calculated ab initio without the need to 
assume an s-like conduction band wavefunction (i.e., no assumption of a parabolic band). 

The rates of the elastic scattering mechanisms, z/ e /, are calculated from the electron group 
velocities, v, and density of states, D s ; thus, the mobility may be calculated directly from 
the electronic band structure. The original form of these equations from k-p theory, and the 
modified equations that we propose, are listed in Table |TJ Note that in every equation, ; 
which, in semi-empirical models, is the group velocity fitted to experiment by the band gap 
and effective mass of the semiconductor (included in d(k), see Table]]]), has been replaced 
by its ab initio counterpart, or v (k), which is calculated directly from the band structure. 

As an example, the DFT-calculated density of states (DOS) can be plugged into Equa¬ 
tion A7 to obtain the inverse charge screening length, /3, in ionized impurity scattering. 
Furthermore, the numerator and denominator of the integrand in Equation [5] both contain 
the density of states of a free electron gas, . Since this can also be calculated ab ini¬ 
tio for the specific system of interest, Ds can instead be substituted in the equation for 
calculating the mobility and reformulated in terms of the energies, e\ 


1 f v (e) D s (e) g (e) de 
3E f D s (e) f (e) de 


where, again, v (k) is the electron group velocity and g is the perturbation to the electron 
distribution, which is calculated iteratively using Equation [2] and can be expressed both as 
a function of k or £ (k) (i.e., the band structure). 

Once the mobilities of the electrons and holes are known, the electrical conductivity can 
be readily calculated: 

a = nep e + pep h (7) 


where n and p are the concentration of electrons and holes, respectively, e is the absolute 
value of the charge of an electron and p e and ph are the mobility of electrons and holes 
respectively. 

It should be noted that there are fundamental differences between the model that we 
have presented here and those relying on the relaxation time approximation (RTA), and 





TABLE I. The original equation^ 3 ®, based on k • p theory for elastic scattering rates and overall 
drift mobility, and proposed modifications, based on ab initio parameters, introduced in this work. 


k • p theory with empirical parameters 


ab initio 


vu (k) = [D (k) In (l + - B (A:)] m (k) = 


e 4 V 


£ 2 = wr/(;) /(i-/)d* 


■V W = [3 - 6c 2 (k) + 4c‘ (*)] 

|]c 2 (*) = l-^, a 2 (i) = H-3|f (^-1) 

( k) = [3 - 8c 2 (k) + 6c* (*)] 

fu\ _ N dis e 4 md(k) _1_ 1 _ i , m/m *-1 

Vdis W — fi a e 2 c 2 ( „ t2 \ 3 /2_ ’ d(fc) — Q 


ta 2 2 / ,\3/2 

n tnC ’ ' _4fc 2 \ , 


( 1+ ^7 


_ K f k 3 (g(k)/Ed{k))dk 
Coverall ~ 3m j k 2 fd k 

g(k) = f (k) - f 0 ( k ) 


SttEqH 2 k 2 v(k) 


D 


(k)\n{l + f-)-B(k) 


^pe 


P 2 = ^ BT jD s (s)f(l-f)de 
c (k) : obtained directly from wavefunctions 

^ ( fe ) = t 3 - 8c2 ( fe )+ 6c4 ( fc )] 


u dis (k) — h 2 e 2*l,a 


n^v(k) 




_ l / «(e)-Ds(e)fl(e)«fc 

Coverall ~ 3E j Ds{e)f(e)de 

g(e) = f (e) - /o (e) 


a The subscripts stand for: ii (ionized impurity), pe (piezoelectric acoustic phonon), de (deformation), and 

dis (charged dislocation scattering). The parameters are: m (electron mass), m* (effective mass), eo 

(low-frequency dielectric constant), e g (band gap), v ( k ) (electron group velocity), D$ (e) (ab initio 

calculated density of states), c(k) (contribution of p-type orbitals to the conduction band), /3 (inverse 

ionized impurity charge screening length), E e (deformation potential), Cj (lattice constant), E (small 

electric field), and c e i (spherically averaged elastic constant). B ( k) and D ( k ) are just collection of the 

parameters: c, k and /3. Their purpose is to simplify the equation 3 . 
b The c(k) parameter is the contribution of the p orbital to the wavefunction of the band. In the k • p 

formulation, it has a closed-form expression that includes the band gap and experimental effective mass. 

In the ab initio formulation, this wavefunction admixture can be calculated by projecting the 

wavefunctions onto spherical harmonics that are nonzero within the sphere around each ion; this 

procedure is already implemented in the Vienna ab initio Simulation Package (VASP^P 1 -®!!. 


particularly, BTE-cRTA. Rather than simplification of the collision term in the BTE (Equa- 
A2) through the RTA (Equation |A3[ ), nlly involve this term by considering both 


tion 


elastic and inelastic scattering mechanisms. It is noteworthy that the BTE-cRTA formula¬ 
tion only implicitly takes into account clastic and inelastic scattering mechanisms, by fitting 
the overall relaxation time to experimental data with no explicit consideration of changes 
in electron distribution from each type of scattering mechanism. Furthermore, unlike the 
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semi-empirical models that were described above, we use ab initio parameters; thus, higher 
predictability and little to no dependence on experimental data is achieved. 

B. ab initio Parameters 

The main input that is needed for the transport model is the crystal structure of the 
semiconductor material, from which ab initio parameters, such as the (optimized) lattice 
constant, PO phonon frequency, dielectric and piezoelectric constants, deformation potential 
and effective mass, can be computed. 

We also need to know the Fermi level to compute scattering rates in Table [T] In order 
to obtain the Fermi level, we calculate the carrier concentration and match it to the given 
concentration (input), n, according to Equation [8] 

1 r- t-OO 

n = y J £ 9 ( £ ) f ( £ ) de ( 8 ) 

Since both of the III-V semiconductors considered here are n-type, the concentration of hole 
carriers is negligible. The concentration of ionized impurities, N u (see Table [T]), is the sum 
of the concentration of all ionized centers regardless of the sign of their charge, since they 
are scatterer centers in both cased^ 

Ab.- 

N H = N a + N d + ^ (9) 

Cl 

where Nd and N a are concentration of donors and acceptors, respectively. N n can then 
be calculated at a given electron concentration, n, by iteratively solving the charge balance 
equation 8 : 

Ab- 

n + N a = N d + — (10) 

Cl 

where the density of dislocations, N^ ls , is only relevant for InN and is considered to be zero 
for GaAs. In both GaAs and InN, temperatures lower than 20 K need not be considered due 
to the deionization of shallow donors at lower temperatures, as observed experimentally 26 . In 
the case of InN, electronic scattering from existing linear charged dislocations thus becomes 
important. The density of the dislocations, N^is, can be determined from TEM images, in 
the units of (cm -2 ). We can thus obtain the overall density in bulk, by assuming that these 
linear dislocations are uniformly developed along the c-axis. This is reflected in dividing 
the dislocation density by the lattice constant, q, in Equation [TOj By doing that, we are 
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assuming that there is one unit of positive charge (donor) per unit cell. For InN samples, 
according to Miller et alP, one can assume full ionization of the donors, and therefore, 
a compensation level of one (i.e., Nl >+ Na = 1 or Njj » AG). Also, the assumption of 
donor or acceptor charged dislocations yields similar result^, therefore, we assume donor 
dislocations dominate here. It should be noted that we compare the calculated N dls with the 
corresponding experimental data if available; otherwise, the limit for electronic properties 
at different values of N dts can be calculated without the need for experimental data. 

On the other hand, in a pure, epitaxially-grown, high-mobility GaAs sample with an elec¬ 
tron concentration of n = 3.0 x 10 13 , no dislocations exist (i.e. AGs = 0). The concentrations 
of donors and acceptors have been separately reported®^, so this provides validation of the 
accuracy of our model, without needing to solve for N tl . However, in the general case where 
the electron concentration is unknown, we can plot the mobility and Seebeck coefficient at 
different compensation ratios to define the limit of the transport properties, as shown in 
Figure [5j Therefore, it is important to note that only when comparing with experimental 
mobilities/Seebeck coefficients do we use experimentally measured electron concentrations; 
otherwise, we may calculate ab initio mobility or Seebeck coefficient, for example, at various 
electron concentrations, without any reliance on experimental data (e.g., as shown in Figure 

§■ 

We use Brooks-Herring theory for singly-charged ionized impurity scattering 2 ^, as shown 
in Table |l| This is supported by the fact that in GaAs, oxygen impurities, O^ 3 , have been 
confirmed to be dominant and singly charged^, while in InN, nitrogen (donor) vacancies, 
Vtf, are dominant and singly charged 28 . It should be noted that the Brooks-Herring formu¬ 
lation is more accurate at low carrier concentrations, since at high concentrations, despite 
the inherent assumption of the theory, not all electrons are screened by the charge of an 
ionized center. More information on the Brooks-Herring ionized impurity model is available 
in Appendix |A 2| 

In order to calculate the low- and high-frequency dielectric constants, we use density 
functional perturbation theory (DFPT), as implemented in VASP, to determine Born effec¬ 
tive charges, dielectric and piezoelectric tensors, including local held effects in DFT, as well 
as the force-constant matrices and internal strain tensors. We then subtract the ionic con¬ 
tribution to the static dielectric tensor to obtain the high-frequency dielectric constant 29 30 . 
Furthermore, the inelastic scattering effect is strongly dependent on the longitudinal polar 
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optical phonon frequencies, u po . These frequencies can be calculated using the Phonopy 
code®, where we identify the highest energy peak in the optical phonon density of states. 
It should be noted that at and around the T point, the phonon frequency is almost constant 
(Figure [8j) . 

To calculate ab initio the deformation potential, Ed, we strain the system and calculate 
the energy of the conduction band of InN and GaAs unit cells at different volumes. Then, 
we approximate the deformation potential using the following equation: 


E n = -V 


'dE, 


CBM 


dV 


( 11 ) 


at V=V 0 


where V is the volume, Eqbm is the energy of the CBM and Vo is the volume of the 
relaxed structure (i.e., zero pressure) 32 33 . It should be noted that since the absolute value 
of Ecbm is a function of the volume itself, we use the difference between the energy of 
the first conduction band and the first valence (core) band. Furthermore, the elastic and 
piezoelectric constants have been already calculated ab initio for GaAs and InN, and are 
available in the literature. For GaAs, we use the values calculated by Beya-Wakata et ah®, 
and for InN we use the values calculated by Sarasamak et ah®, to obtain the piezoelectric 
coefficient and elastic constant used in the equations for piezoelectric scattering in Table [T| 
As a comparison, the electrical conductivity and Seebeck coefficient are also computed 
using the widely-used BTE-cRTA formulation. We choose the BoltzTraP package®, which 
uses Fourier interpolation of the calculated bands, and differentiate the band energies to 
find the group velocities of the electrons. Other than the need to fit the relaxation time to 
experimental measurements of the conductivity, the BoltzTraP/BTE-cRTA implementation 
represents an otherwise parameter-free model that can be adapted to different semiconductor 
materials. 


III. COMPUTATIONAL METHODOLOGY 

For each semiconductor material, the geometry of the unit cell is optimized, and the 
density of states and band structure are calculated. In the case of zinc blende GaAs and 
wurtzite InN, the unit cells are optimized using Kohn-Sham density functional theory (KS- 

dft ^3BI3?| 

as implemented in VASP. The generalized gradient approximation of Perdew, 
Burke, and Ernzerhof (GGA-PBE) 38 39 is used to express the exchange-correlation potential, 
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TABLE II. Structure of GaAs and InN calculated with DFT, using the GGA-PBE exchange- 


correlation functional. Changes in the lattice constants compared to experimental valued 2 43 upon 
optimization are reported below. 


Compound 

Space Group 

l a (A) 

% change in a 

l c l (A) 

% change in c| 

GaAs 

F-43m 

5.75 

2.17% 

- 

- 

InN 

P 631 UC 

3.533 

0.56% 

5.693 

0.8% 


and Projector Augmented Wave (PAW) potential^ 40 41 are used to represent the valence 
wavefunctions. Information regarding the structure of these two systems and their changes 
upon geometry optimization have been summarized in Table [TT[ The initial structures are 
obtained from the literature 42 43 . 

We then compute the electronic band structure of these materials. The energy cutoff 
for the plane wave basis set is set to 500 eV. The band structure is computed in line 
mode along seven high-symmetry fc-points in the IBZ, with 20 Appoints between each pair of 
high-symmetry points. The self-consistent density of states (DOS) calculation is performed 
using a 20 x 20 x 20 fc-point mesh, for both GaAs and InN. The non-self consistent energy 
calculations are performed in a special k-point mesh around the T point, at which the 
conduction band minimum (CBM) occurs in both direct band gap GaAs and InN. This 
fc-point mesh contains a total of 10,234 points in the Irreducible Brillouin Zone (IBZ), with 
mesh spacing of 0.001, 0.01, or 0.1 fractional units, to completely account for band anisotropy 
while remaining dense enough around the T point to obtain accurate group velocity and 
effective mass values. To determine the effect of presumably more accurate band structure 
calculations on the band curvature, effective mass, and group velocity, we have also employed 
the GW method. Only 941 fc-points in the IBZ have been used for GW calculations, since 
this method is more computationally demanding. Using more fc-points does not change the 
calculated effective mass. The GW 0 band structure calculations are performed using the 
maximally-localized Wannier functions (MLWFs) interpolation, as implemented in VASP 
and the Wannier code 44 . It is very important to also show the feasibility of the ab initio 
model with band structures calculated using DFT+U, since the GW and hybrid functional 
(e.g., HSE 45 ' 4 ') methods are computationally demanding for complex materials with larger 
unit cells than GaAs. On the other hand, for InN, all methods and functionals attempted, 
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including LDA, GGA, HSE and GWO, resulted in a band structure with zero band gap and 
a falsely predicted p-like conduction band. Only GGA+U®!, with U values obtained from 
the literature 49 50 (U^ = 6 for In and U p = 1.5 for N), produced a correct band structure and 
with a more s-like conduction band particularly around T point, which is consistent with the 
self-interaction corrected band structure reported by Furthmullcr et al. 51 . In the process of 
choosing the U value for GGA+U band structure calculations on GaAs, however, the values 
(U = 8 eV for both d orbitals of Ga and As) recommended by Persson and Mirbt 5 -^, with 
an emphasis on correctly obtaining the band gap and effective mass values, result in GaAs 
falsely becoming an indirect semiconductor, with the conduction band minima located at 
the L and A" /c-points rather than the T point 52 . Therefore, we have also employed effective 
U values of 7 eV (Ga) and 6 eV (As), for which a direct band structure is obtained. We have 
calculated mobilities obtained from both of these band structures and compared them with 
the ones obtained by the GW band structure. In order to calculate the group velocities, 
v ( k ), and the overall average effective mass, we have fitted a sixth degree polynomial to the 
calculated conduction band (i.e., average energy as a function of distance from T point or 
e(k)) with R 2 > 0.99: 

(12) 


m 


* 


{ 1 d 2 e 
yft 2 dk 2 


-i 


at k =0 


(13) 


It should be noted that we do not use the value of effective mass in the proposed carrier 
transport model. Rather, we calculate it solely to compare with experiment and evaluate 
the effect of the shape of the conduction band (i.e., group velocities) calculated by various 
methods, such as GGA, GGA+U, and GW. Fitting polynomials to the numerically calcu¬ 
lated conduction band and density of states results in smooth plots of mobility and Seebeck 
coefficient, as presented here, while preserving the values that are calculated ab initio with 
R 2 > 0.99 in all segments fitted. We fit these polynomials at different segments of the band 
structure and carefully choose only the ones that result in the maximum R 2 and minimum 
discontinuity where the polynomials meet. This results in very smooth calculated group 
velocities, and, subsequently, other transport properties. 
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(a) GaAs 

FIG. 1. Band structure of cubic GaAs and wurtzite InN, normalized so that the Fermi level is set 
to zero at the conduction band minimum. 


IV. RESULTS AND DISCUSSIONS 


A. ab initio Calculated Parameter Inputs to the Transport Model 

The computed band structures of GaAs and InN are shown in Figure [lj We have cal¬ 
culated a G1U0 band structure, which starts from the wavefunctions previously computed 


using the GGA-PBE functional, as shown in Figure la 


The band structures used in previous semi-empirical modelk^® express the energy of the 
conduction band as a function of the distance from the Y point. Instead, we calculate the 
ab initio band structure in a three-dimensional grid around the CBM, and then average the 
energy values of the ^-points that share the same distance from the Y point (Figure [2]). For 
both GaAs and InN, the ab initio and k • p band structures agree well at small /e-points; 
however, they diverge at larger fc-points. This directly impacts the group velocity of the 
electrons and, ultimately, the transport properties - particularly at higher temperatures 
where higher energy electrons have nonzero occupation. 

We have also calculated a GGA+U 48 band structure, with U values taken from the pub¬ 
lished literature 49 -®, as shown in Figure [lj For InN, GGA+U correctly yields an s-like 
conduction band and a band gap of 0.5 eV, which is comparable to the self-interaction cor¬ 
rected band gap of 0.58 eV reported by Furthmuller et al. 51 and the experimental values of 
0.675-0.7 e\® 57 ® (Table |hT]). We include DFT+U calculations only to show the feasibility 
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TABLE III. Inputs to the transport model, as calculated ab initio compared to experimentally 
measured values. The bolded numbers are used in our transport property calculations; note that 
not all appear in the final expressions listed in Table [I] 


Parameter 

GaAs 

InN 

ab initio 

Exp. 

ab initio 

Exp. 

ci (nm) 

0.562 

0.575® 

0.565 

0.569® 

uipo (THz) 

8.16 

8.73® 

17.83® 

17.65 

eo 

12.18 

12.91® 

11.42 

10.3® 

^OO 

10.32 

10.91® 

6.24 

6.7® 

E d (eV) 

6.04 

8.6® 

4.46 

3.6® 

>k 

m 

0.053-0.066 

0.0636-0.082 EMU 

0.062, 0.071 (GW) 

0.05-0.08 ES® 

s g (eV) 

0.96, 1.19 (GW) 

1.424® 

0.50 

0.675-0.7 E 57 61 


a The parameters are: c; (lattice constant), w po (PO phonon frequency), eo (low-frequency dielectric 
constant), (high-frequency dielectric constant), Ed (deformation potential), m* (effective mass), e g 
(band gap). 

b The GaAs effective masses are calculated as 0.053 (GGA+U, this work), 0.066 (GGA+U, with published 
iM), and 0.063 (GW0). 


of these less-expensive methods, in the case of more complex semiconductor materials for 
which a GW calculations is too expensive. Also, DFT usually suffers from vastly underesti¬ 
mating the effective mas£p-®, and the introduction of the fitting parameter U may reduce the 
predictability of the ab initio model as a whole. Therefore, we stress that all reported trans¬ 
port properties are calculated here using the parameter-free GW band structures, unless 
otherwise stated. 

Although we do not directly use the value of the electron effective mass in the transport 
property expressions, we see that the calculated effective mass of 0.062 for InN is consistent 
with the previously calculated effective mass (0.066) using an empirical pseudopotential® 2 , 
and well within the range (0.05-0.08) measured experimentally 55 60 . 

We also show the calculated phonon band structure and density of states of these two 
compounds in Figure [8j For GaAs, the calculated PO-phonon frequency of 8.16 THz is 


shown in Figure 8a For InN, the calculated optical phonon frequency of 17.83 THz is close 
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k (nm" ) 

(b) InN 

FIG. 2. The conduction bands expressed in terms of the average energy as a function of distance 
from the CBM (i.e., center of Brilloun zone, or T point), as calculated from semi-empirical expres¬ 
sions (in k ■ p formulation) versus ab initio. The difference at higher k values has a significant 
impact on transport properties, especially at high temperatures. The values of U for the d orbitals 
of Gallium and Arsenic, respectively, are in parentheses, while those for InN are taken from the 
published literature 49 50 . 


to the 17.65 THz value reported by Bungaro et al. 55 63 . We have listed all the parameters 
that are used in our transport model in Table [TTTj We have calculated all of these parameters, 


as bolded in Table III, ab initio to demonstrate the feasibility of a fully predictive model for 
transport properties. The only exceptions are the elastic and piezoelectric constants, which 
are necessary to calculate the piezoelectric coefficient, P, in Table [[} As described earlier, 
we have instead used the previously calculated values from published DFT studies for these 
constant ^ 34 1 35 ! 


B. Model Validation on GaAs 

In order to evaluate the accuracy of our model, we first calculate the mobility of three 
experimentally synthesized and characterized GaAs samples, as described by Stillman et 
al. 26 . We also perform this analysis over a wide temperature range for high purity GaAs 


samples with very low electron concentrations, as labeled as ’’pure” in Table IV 
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TABLE IV. Carrier concentrations of various experimentally fabricated and characterized GaAs 
samples. For the ’’pure” sample, data is available roughly between 5-1000 K. For the real samples, 
mobility data is also tabulated at different temperatures. 


Sample 

Concentration, n (cm 3 ) 

Donor, Np 

Acceptor, Na 

Reference 

pure 

3 x 10 13 

5.2 x 10 13 

2.2 x 10 13 

El 

a 

2.7 x 10 13 

4.8 x 10 13 

2.1 x 10 13 

ED 

c 

7.7 x 10 14 

1.1 x 10 15 

3.3 x 10 14 

ED 

e 

3.1 X 10 15 

4.7 x 10 15 

1.6 x 10 15 

ED 


As shown in Figure [3a| the most accurate GW band structure results in the best agree¬ 
ment with experimental data. The DFT+U band structure, however, does provide us with 
limits of the mobility over different temperatures. When calculating the mobility and See- 
beck coefficient, we calculate the Fermi level by first calculating the electron concentration 
through Equation [8j and then matching it to a given concentration. The calculated prop¬ 
erties are very sensitive to the calculated Fermi level. Therefore, for comparison, we have 
included the results using both the ab initio DOS used in Equation [8j and the free electron 
DOS. As shown in Figure [3a] the ab initio model for DOS performs better for lower electron 
concentrations and lower temperatures, while the free electron DOS is more suitable for 
higher temperatures, and, particularly, at higher electron concentrations. We acknowledge 


that because of the log scale in Figure 3a, seeing the quantitative agreement is difficult. 
Therefore, we report the calculated relative error compared to the experiment for the best 
cases for each sample - from the ab initio DOS for sample a and from the free electron 
DOS for samples c and e. The minimum, maximum and the relative error in calculating 
the mobility of sample a are 2.25% (at 195 K), 29.42% (at 29 K), and 13.33%, respectively. 
These numbers are 1.02% (at 167K), 15.01% (at 49K), and 7.97% for sample c and 0.22% (at 
195K), 7.90% (at 40K), and 4.04% for sample e. Overall, the agreement is poorer at higher 
electron concentrations and lower temperatures; this is attributed to the inaccuracy of the 
Brooks-Herring ionized impurity scattering model at high electron concentrations, as briefly 
described in Section [TT| Furthermore, the model has also been validated with the data on 
crystalline samples with very high purity. The calculated electron mobilities, assuming the 
limit that only one scattering mechanism exists at a time, along with the overall mobility, 
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temperature and electron concentration to 


experimental data. The values of U used for the 
d-orbitals of Gallium and Arsenic, respectively, are 
in parentheses. 



Temperature (K) 

(b) Limitation of ’’pure” GaAs mobility from each 
scattering mechanism. 


FIG. 3. The calculated and experimental 64 mobility data for GaAs at various electron concen¬ 
trations and temperatures. More details on the experimental data, including donor and acceptor 


concentrations, are available in Table IV 


are shown in Figure [3b] The reasonable agreement between the calculated and experimental 
mobilities provides independent validation of the transport model. The minimum, maximum 
and average relative error of calculated mobility are 0.46% (at 394K), 23.55% (at 175K) and 
9.53% respectively for temperatures above 20 K. The mobility is mainly limited by ionized 
impurity scattering at low temperatures, piezoelectric scattering at intermediate tempera¬ 
tures, and polar optical phonon scattering at higher temperatures (> 60 K); all of these are 
consistent with the previous results shown by semi-empirical models™ ® yet no experimental 
parameter has been used here in predicting the correct changes with the temperature and 
the carrier concentration. 

Once we have the calculated mobility, at a given electron concentration, we can calculate 
the electrical conductivity of GaAs by Equation [7| For now, we assume that the carrier con¬ 
centration remains constant with temperature over the range of interest. We then compare 
to the experimental conductivity and those values calculated using the BTE-cRTA frame- 
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Temperature (K) 

FIG. 4. Electrical conductivity of GaAs calculated using the ab initio transport model and the 
BTE-cRTA framework, and compared to experimental data. The Fermi level is calculated by 
matching the calculated carrier concentration to n = 3 x 10 13 . This has been done either at the 
mentioned temperature and kept constant over the whole temperature range, or in the case of 
’’matched Fermi”, at each temperature, the Fermi level is adjusted to the given n. The relaxation 
time, r, is determined by fitting the calculated conductivity to the corresponding experimental 
value at 300 K. The calculated value of 4.5 x 10~ 23 s for t is unreasonably low but it has been 
included for the sake of comparison of all models. 

work, under the scenarios listed in Figure [4j As shown, not only does BTE-cRTA fail to 
correctly predict the trend for conductivity with temperature, but also quantitatively differs 
from the experimental values. 

Finally, we calculate the Seebeck coefficients of the GaAs samples (assumed to be at 300 
K), and compare them to the values reported previously by Rode and Knight 4 (Figure [5]). 
Since the data are for various samples with different electron concentration and compensation 
ratios, we choose various values of Na/n = (Np + NjCj/n. As shown, a range of Seebeck 
coefficients are calculated at each electron concentration, which includes the experimentally 
measured points. It should be noted that not knowing beforehand the compensation and 
concentration of donors and acceptors, as well as their charge states, limits the overall 
predictability of our model. However, even given these limitations, the close fit between 
ab initio and experimental properties provides independent validation of the viability of 
our model. For further evaluation, we have calculated the Seebeck coefficient, assuming 
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FIG. 5. Calculated and experimental^ GaAs Seebeck coefficient, at different ratios 
ized impurity concentration, Nu, to the electron concentration, n. For the Pisarenko 


Equation 14, we have used values of r = — \ for acoustic phonons, and m* = 0.063. 


of the ion- 
plot, using 


Pisarenko behavior and compared it to our model in Figure [5] We use Equation |I4| with two 
fitting parameters: effective mass, m* and r. It should be noted that in the case where the 
best agreement with experiment, through Pisarenko behavior, is only achievable by choosing 
either m* = 0.11 or r = 0.35, both of these values are far from experimental measurements 
and thus lack physical meaning. 




VB 


5 , 2 (27 Tm*kBT) 3 ^ 2 

- + r + In ;Ai 


(14) 


C. Model Validation on InN 

In order to further evaluate the accuracy of our model and its applicability to more com¬ 
plicated semiconductors, we also calculate the mobility and Seebeck coefficient (Figure [6]) 
of three experimentally synthesized and characterized InN samples by Miller et al.P These 
calculations are more challenging due to the reported presence of linear charged dislocations 
in the crystal structure 8 due to the processing conditions employed. For each sample 
at a given carrier concentration, as shown in Table [V], we change the concentration of dis¬ 
locations, Ndis, until the calculated mobility values match the experimental measurements. 
The fitted N dis (Table [V]) is within the range of measured concentrations from transmission 
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FIG. 6. Calculated and experimental transport properties of the three InN samples listed in Table 
[v| The dashed line is calculated by the semi-empirical model used by Miller et alP, while the solid 
line is calculated by the proposed ab initio transport model. 




TABLE V. Measured® and calculated InN dislocation density, corresponding to the mobility and 


Seebeck coefficient reported in Figures 6a and 6b 


Sample 


N d is (cm 2 ) 


Experimental 

Semi-empirical® 

This work 

A 

« 1 x 10 11 

1.5 x 10 11 

8.20 x 10 10 

B 

2 - 5 x 10 10 

1.5 x 10 10 

1.18 x 10 10 

C 

« 1 x 10 9 - 5 x 10 10 

4.1 x 10 9 

3.47 x 10 9 


electron microscopy analysis (TEM) 8 !, which confirms that the limiting mechanism is indeed 
scattering from dislocation lines. 


As shown in Figure 6a and 6b, while there is an excellent agreement between the cal¬ 
culated and experimental mobility, the calculated Seebeck coefficients for samples B and C 
exhibit more pronounced changes with temperature than the experimental Seebeck coeffi¬ 
cients. The mobility of the samples is found to be limited by charged dislocations, particu¬ 
larly at low temperatures. The next limiting mechanism is polar optical phonon scattering, 
which is more important at higher temperatures while ionized impurity scattering is more 
important at lower temperatures. This can be seen in Figure [7J which shows the mobility 
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FIG. 7. Calculated and experimental^- values for InN mobility at n = 9 x 10 17 cm 3 (sample B in 
Table 0. Each line represents the mobility if limited only by the corresponding mechanism. 


of sample B if it were limited by each type of scattering mechanism, as well as the overall 
mobility. These findings are in agreement with the semi-empirical transport model 8 } except 
that all parameters are obtained from ab initio calculations that require knowledge only of 
the crystal structure of the material. Comparing the transport properties calculated from 
using model with those calculated using semi-empirical models (including experimentally 
measured band gap and effective mass (See Table III under ”Exp.”) in Figure [6] shows that 
although quantitative agreement with experiment is slightly better with the semi-empirical 
model, Seebeck coefficient calculations on samples B and C, and the mobility of the sample 
at high temperature, show much better accuracy with the ab initio model presented here. 


Finally, we should once again acknowledge the assumptions and limitations of the current 
model when applied to the other types of semiconductors. Most importantly, the formula¬ 
tion presented in this work is for low-field transport (particularly drift mobility and Seebeck 
coefficient), in which the changes to the electron distribution are merely a linear pertur¬ 
bation to the equilibrium Fermi-Dirac distribution; thus, the applicability of the current 
model for high-field transport or heavily doped and polar semiconductors where the linear 
BTE formulation fails 68 , is very limited. Furthermore, we have averaged the energy around 
CBM and expressed the energy values in the band structure as a function of the absolute 
value of k, or simply, the distance from T point in the reciprocal space. Therefore, the 
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reported mobility values are averaged and the effect of band structure anisotropy is not fully 
captured. It is possible, however, to include the band structure of the material only in the 
specific orientation of interest to account for anisotropy. Currently, the model is limited 
to a single conduction band. Although the single band ab initio model can be used for 
prediction of many direct band gap semiconductors, it will only result in an overestimation 
of transport properties of semiconductors with more complex band structure. This is due to 
the fact that currently, interband scatterings between several bands that are participating 
in transport are neglected. In future, we will solve coupled-BTE and take into account two 
and more participating bands which enables calculation of both electron and hole mobilities 
in more materials. Finally, although the usage of the Hubbard U parameter in the band 
structure calculation might limit the predictability of the model in calculating overall trans¬ 
port properties, this can be properly addressed by using more accurate methods of band 
structure calculations as reported here. We include DFT+U calculations here only to show 
the feasibility of working with the model when GW or other less commonly used methods 
are not technically or otherwise feasible. 


V. CONCLUSIONS 

We have presented an ab initio transport model for calculating the electrical mobility 
and Seebeck coefficient of n-type semiconductors. By using the inputs from density func¬ 
tional theory calculations, and considering all relevant physical phenomena (i.e., elastic and 
inelastic scattering mechanisms), we have successfully calculated highly-accurate transport 
properties of GaAs and InN over various ranges of temperature and carrier concentration. 
Our model provides both qualitative and quantitative improvements in accuracy compared 
to the widely-used semi-empirical and constant relaxation time approximation model solu¬ 
tions to the Boltzmann transport equation. Future work will focus on extending this model 
to p-type semiconductors. 
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Appendix A: Boltzmann Transport Equation 


The Boltzmann transport equation (BTE) describes the non-equilibrium behavior of 
charge carriers (e.g., electrons or holes) by statistically averaging over all possible quan¬ 
tum states. For the electron distribution, /, this is represented by the BTE: 


df(k,r,t) (df(k,T,ty 


f] k 

- ^ ■ Vfc/ (k, T, t) - V (k) • V r / (k, T, t) (Al) 


dt \ dt 

\ /a 

where / is a function of state k, temperature T, and time t, and v (k) are the electron group 
velocities. The three terms on the right-hand side of Equation |A1| refer, respectively, to 


the temporal rate of change of / due to all scattering processes, rate of change of / due to 
external forces, and diffusion from the carrier density gradient. 

If the external forces consist only of a low electric field, E. and no magnetic field, B, such 
that A = Ay then the low-field BTE becomes: 


+ v (k) . Vr/ (k , T) + f . vu (k, T) = MM 


(A2) 


1. Constant Relaxation Time Approximation 

Furthermore, / can be described as a first-order (linear) perturbation, g(k), from the 
(equilibrium) Fermi-Dirac distribution, f 0 , due to scattering: 

/<9/(k,T,f)\ / (k) - f Q (k) g(k) 


dt 


fo [ £ ( k )] = - 


(A3) 


(A4) 


e [e( k )-e F ]/fc s T _|_ 

where the dependence of £ on k is given by the electronic band structure, and the various 
scattering terms and time dependence are lumped into the electronic relaxation time, r. 

If r is a constant, then this major simplification results in the BTE-cRTA. This assump¬ 
tion simplifies the theory to an extent that closed form expressions for conductivity and 
Seebeck coefficient can be obtained 10 ! In this approach, the details of all elastic and inelas¬ 
tic scattering mechanisms are lumped into the relaxation time constant, r. While popular, 
this approach suffers from the following disadvantages: 1. r is obtained by fitting to the 
experimental data for the conductivity of the material, which limits the predictability of the 
model, and 2. Due to oversimplification of the transport mechanism, it may result in incor¬ 
rect values and even incorrect trends with temperature or carrier concentration, as illustrated 
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in Figure [4j Therefore, by explicitly including all possible electronic scattering mechanisms, 
one can determine which mechanisms are physically relevant for a given semiconductor. 

2. Explicit Solution of Linear BTE 

To go beyond the relaxation time approximation, both elastic scattering mechanisms, for 
which the kinetic energy of the electrons remains constant, and inelastic scattering mecha¬ 
nisms, for which there is a change in the electron distribution, should be taken into account. 
If the system is governed only by elastic scattering mechanisms, the relaxation time, r, is 
equal to the inverse of the overall elastic scattering rates, which is the sum of all individual 
rates. Evidently, r is not constant but does depend on energy; however, it does not necessar¬ 
ily follow a power law dependence (e.g., in InN®). However inelastic scattering mechanisms 
also limit the mobility, and therefore, the conductivity, of the semiconductor; as an exam¬ 
ple, polar optical (PO) phonon scattering is the main electron-phonon interaction that limits 
mobility at high temperatures in GaAs. Thus, we need to first calculate the perturbation, 
g, to the electron distribution due to elastic and inelastic scattering mechanisms, and then 
integrate g over all states to obtain the mobility. Details on this approach are given below. 

The most relevant elastic scattering mechanism for compound semiconductors is expected 
to be ionized impurity scattering at low temperatures. Ionized impurity scattering occurs 
when a charged center is introduced inside the bulk material. As a result of Coulombic 
interactions between the electron and ion, electrons scatter to different states (i.e., become 
distracted). The ionized impurity scattering rate, Ua (i.e., a component of the overall v), 
may be expressed using Brooks-Herring theory 25 : 

e A N 

8TTelh 2 k 2 v 

where the charge screening potential, 0, is obtained by solving Poisson’s equation: 

0= 4 ^ eXp(_/3r) (A6) 

and inverse screening length, (3, is given by: 

0* = jjA / D s (e )/(!-/) ds (AT) 

where / is the electron distribution and eo is the low-frequency dielectric constant. Details 
on the a, D and B parameters are given in the literature 3 . 


Dl + 


(A5) 
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At high temperatures, after an inelastic (e.g., electron phonon) scattering event, where 
the electron scatters from momentum state k to k', the energy of an electron changes, 
and hence, the electron distribution also changes. (Note that the distribution may also be 
perturbed by external forces, such as an electric field or temperature gradient.) Thus, / 
becomes a function of k, so it must be mapped via the electronic band structure, e (. k ). This 


effect can be shown as a deviation from Fermi-Dirac behavior (Equation A3). After some 
mathematical manipulation, for which details can be found in the literature^, the BTE can 
be reformulated as: 

Si (s') - - (10 - ¥ u 


9 = 


(AS) 

(A9) 

(A10) 


^el 

Si (sO = f dk'Xg ( k ') [s ine i (. k ', k) [1 - f (k)] + s ine i (. k , k’) f (k)] 

S 0 = I dk' [Si ne i {k, k') [1 - / ( k')] + s ine i (k', k ) / ( k ')] 

Detailed integrated expressions for the scattering in, S), and scattering out, S 0 , terms are 
available in the literature 3 . The reformulated BTE can then be solved iteratively, using 
Rode’s method® 8 AM3 since g. (r/) and / themselves are functions of g. First, the Fermi- 
Dirac distribution can be plugged into the right-hand side of Equation |A8| to obtain the 
first guess, r/i, which in turn is used to obtain a new electron distribution to solve for the 
next guess, < 72 ; this process continues until g converges to a unique value. Typically, five 
iterations are required for the perturbation to converge for polar optical phonon scattering 
in GaAs and InN. More details on Equations A8 A9 are available in the literature 3 . 


3. Phonon Dispersion 

Polar optical phonon scattering originates from interactions between electrons and high- 
frequency optical phonons. They provide the dominant inelastic electron scattering mecha¬ 
nism near (and above) room temperature in compound semiconductors. This is attributed 
to the high energies of optical phonons being comparable to at high temperatures. The 
scattering rates themselves are strongly dependent on the polar optical phonon frequencies. 
l o po . These frequencies can be calculated using the Phonopy code 31 which solves for dy¬ 
namical matrix from the force constants calculated using density functional perturbation 
theory (DFPT), as implemented in VASP. The phonon band structures for GaAs and InN 
are shown in Figure |8j 










(a) GaAs (b) InN 

FIG. 8. Phonon band structures of InN and GaAs calculated by Phonopy 31 . 

Appendix B: Sensitivity analysis 

1. Sensitivity to the calculated dielectric constants 

We have performed a sensitivity analysis for the calculated mobility of the GaAs pure 
sample at different dielectric constants. As shown in Figure [9j the result is sensitive to 
dielectric constants at low and high temperatures but much less sensitive at temperatures 
in the 100-200 K range. Inaccurate calculation of dielectric constants by -20% can result 
in errors of up to -41% (at 40 K) in the calculated mobility, compared to experimentally 
measured values; for deviations of +20% in the dielectric constants, the resulting mobility 


can increase by up to +43% (at 5 K). The base values are the ones reported in Table III 


as calculated ab intitio and assuming the relaxed structure. This shows the importance of 
accurate calculation of these constants with at least 5-10% accuracy. 


2. Sensitivity to the lattice constants 

We also applied ±3% strain to the lattice constant of the relaxed GaAs unit cell, and 
recalculated the band structure, DOS and optical phonon frequencies to ascertain the effect 
on the calculated mobility. We assumed that everything else is kept constant according to 
the base values (see Table III). According to Figure [lOj the calculated mobility is extremely 
sensitive to the crystal structure. This is mainly due to the impact that the structure has on 
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Temperature (K) 


FIG. 9. Sensitivity analysis of the mobility of the GaAs pure sample (see Table IV). We changed 
here only the static, e s , and high frequency, £oo, dielectric constant from -20% to +20% from the 


base values reported in Table III The results are sensitive to dielectric constant at low and high 
temperatures. 


the band shape (i.e. group velocity of the electrons) since the mobility at any temperature 
is affected. For example, the GW band structure of -3% strained GaAs gives an effective 
mass of 0.026 while that of +3% strained GaAs gives an effective mass of 0.10. Both of 
these values are well outside of the range of the reported experimental values (0.064-0.082, 


see Table III). Also, these strained structures are extremely unlikely to be relaxed with any 
functional, since their built-in pressure with GGA-PBE functionals are already 10.66 kB and 
-74.77 kB, respectively, while the relaxed structure that we calculated and reported in Table 
[TT| has a built-in pressure of only -0.3 kB. Nevertheless, Figure 10 shows the importance 
of accurate calculation of the crystal structure, and subsequently, the band structure (i.e., 
group velocities). 


30 













Temperature (K) 


FIG. 10. Sensitivity analysis of the mobility of the GaAs pure sample (see Table IV). We changed 


here the crystal structure, and subsequently, the newly calculated optical phonon frequencies. The 


calculated mobility is sensitive to the strain at all temperatures. 
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